The skyrmion lattice phase in three dimensional chiral magnets from Monte Carlo 

simulations 
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Chiral magnets, such as MnSi, display a rich finite temperature phase diagram in an applied 
magnetic field. The most unusual of the phases encountered is the so called A-phase characterized 
by a triangular lattice of skyrmion tubes. Its existence cannot be captured within a mean-field 
treatment of a Landau-Ginzburg functional but thermal fluctuations to Gaussian order are required 
to stabilize ivr*. In this note we go beyond Gaussian order in a fully non-perturbative study of 
a three dimensional lattice spin model using classical Monte Carlo simulations. We demonstrate 
that the A-phase is indeed stabilized by thermal fluctuations and furthermore we reproduce the full 
phase diagram found in experiment. The thermodynamic signatures of the helimagnetic transition 
upon cooling from the paramagnet are qualitatively consistent with experimental findings and lend 
further support to the Brazovskii scenario^ which describes a fluctuation driven first order transition 
due to the abundance of soft modest 



I. INTRODUCTION 



Chiral magnets like MnSi or Fei-^Co^Si have received 
a lot of interest recentlyi^"^, mainly by virtue of them 
showing a thermodynamic phase which is characterized 
by a lattice consisting of tubes of magnetic skyrmions. 
Besides the very existence of this phase there seems to 
be huge potential to use these materials for spintronics 
applications^. 

The lack of inversion symmetry in the crystalline struc- 
ture of these magnets gives rise to weak Dzyaloshinskii- 
Moriya (DM) coupling. The competition of this inter- 
action with the much stronger ferromagnetic exchange 
(FM) results in a twist in the magnetic order, leading to 
helical order. Since the DM coupling is weak compared 
to the FM exchange coupling there are long modulation 
periods of many lattice constants, e. g. in the chiral mag- 
net prototype MnSi the modulation period is about 190 A 
while the lattice constant is only 4.6 A^. The competi- 
tion between these two types of interactions determines 
the length of the magnetic spirals but not their direction. 
Consequently, one expects a large ground state degener- 
acy at zero magnetic field. This degeneracy is, however, 
lifted by weak crystal anisotropies which provide an easy 
axis for the ordering wave vector (e. g. [Ill] in MnSi). As 
a direct consequence the phase with helical order has a 
single ordering wave vector. If additionally a finite mag- 
netic field is applied it becomes energetically favorable to 
have the ordering wave vector point in direction of the 
magnetic field. All spins then point in a plane perpendic- 
ular to the field and the system can gain energy by simply 
tilting all spins continuously out of that plane in direction 
of the field, leading to a spiralling umbrella structure. 
This state is referred to as conical phase. Depending on 
the direction of the field the phase transition between 
these two phases is either first order or a crossover and 
occurs at some field value B c \ where the energy gain from 
tilting all spins towards the field becomes larger than the 
crystal anisotropy energy difference between the two di- 
rections of the ordering wave vector. Figs. fl]a.) and b.) 



schematically show the magnetization structure in these 
two phases. 

In 2009 neutron refraction experiments on MnSi 1 dis- 
covered a new thermodynamic phase at intermediate 
fields and temperatures just below T c ss 30 K. This phase 
is characterized by a periodic arrangement of tubes of 
skyrmions which arrange on a triangular lattice (see Fig. 
[T]c.)). Consequently, this phase is referred to as skyrmion 
lattice phase or A-phase. The skyrmion lattice phase 
can be pictured as a superposition of three helices with 
equal pitch length and relative angle of 120 degrees in 
the plane perpendicular to the magnetic field. While in 
mean-field theories based on a minimal Landau-Ginzburg 
theory for anisotropic non-centrosymmetric magnets the 
skyrmion lattice was argued to be a stable solution, this 
phase does not appear as a stable phase and always is 
slightly higher in its free energy than the conical phase 
for cubic systems^, such as MnSi. While it was argued 
that this phase could still be stabilized by long-ranged in- 
teraction d 10 * 11 * or extra phenomenologica]_ parameters^ 
in the free energy functional, MiihlbaueAl et al. found 
that a very natural alternative mechanism to stabilize 
the skyrmion phase is given by thermal fluctuations to a 
Gaussian level on top of the mean-field theory. 

In order to make this argument stronger it is desir- 
able to use an approach which is not based on Gaussian 
fluctuations, but instead incorporates the thermal fluc- 
tuations in a non-perturbative manner, namely classical 
Monte Carlo (MC ) sim ulations. Simulations for two di- 
mensional system#E31 have been performed before and 
indeed found a stable skyrmion lattice phase. The phase 
diagrams obtained are in excellent agreement with recent 
experiments^ on thin films of Feo.5Coo.5Si even though 
the itinerant character of the underlying electronic sys- 
tem is not taken into account in these studies. A ma- 
jor reason why these studies have not been extended to 
three dimensional systems yet is the high computational 
effort: large system sizes are required to account for the 
long spatial modulations. A further complication stems 
from the fact that one has to be very careful in choosing 
effective parameters of the underlying lattice-spin model. 



In this paper we fill this void and perform a MC study 
for three dimensional chiral magnets. We confirm that 
the effect of thermal fluctuations indeed is what stabilizes 
the skyrmion phaseM Overall, we find excellent agree- 
ment with the experimentally observed phase diagram as 
well as with non-trivial thermodynamic signatures across 
the phase boundary from the paramagnet into one of the 
respective ordered phases. For zero magnetic field the 
transition from the paramagnet into the helimagnet is a 
fluctuation-driven first order transition and can be de- 
scribed in terms of the Brazovskii scenario^-. Some of the 
experimental features of this transition have proven hard 
to capture in purely analytical approaches 3 , however, the 
Monte Carlo approach captures all the qualitative fea- 
tures. 

The organization of the paper is as follows: We start 
with a discussion of the model and the method. Most im- 
portantly, we introduce a minimal lattice model which is 
consistent with the system symmetries and consequently 
the Ginzburg-Landau functional. We furthermore intro- 
duce the Metropolis algorithm together with the algo- 
rithm which is required to overcome the large hystere- 
sis in the underlying system. From there we move to 
the global phase diagram of a chiral magnet in an ap- 
plied magnetic field and compare it to experimental find- 
ings. We close with a comparison of some thermody- 
namic quantities to the experimental findings in zero and 
non-zero field as we go across the thermal transition from 
the high temperature paramagnet towards one of the or- 
dered phases. We find excellent agreement with experi- 
ment which lends further support to the relevance of our 
approach to this problem, despite the rather small lattice 
sizes. 



II. MODEL & METHODS 

A. The lattice Hamiltonian 

Assuming a slow variation of the spin textures one can 
resort to a coarse grained continuum model for the de- 
scription of the magnetic properties of chiral magnets. 
The commonly used one assumes the form (cfP^l) 
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consisting of ferromagnetic exchange J, magnetic field 
B, and a DM interaction K. Above, a is the typical 
distance over which the spin structure can be treated as 
uniformly ordered allowing for the coarse graining pro- 
cedure. This effective model has to be understood in 
connection with the renormalization group meaning that 
terms accounting for the actual microscopic lattice struc- 
ture can be dropped by virtue of them being irrelevant at 



the critical point. They can, however, be important for 
a faithful description deep inside the ordered phase. In- 
stead of the full B20 structure of MnSi one compactifics 
the above continuum theory onto a simple cubic lattice 
(which in principle has inversion symmetry unless explic- 
itly broken as we do below). The construction principle 
is that the effective low-energy theory which can be de- 
rived from the lattice Hamiltonian agrees with the above 
model up to terms which are irrelevant in the renormal- 
ization group sense. 

As a microscopic model we adopt the lattice Hamilto- 
nian presented irPSl, 



H — —J 2_^ S r • (S r +x + S r+ y + S r+ z) — B • 2_^ S r 
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In the following we discuss how to extend this model 
in order to get rid of discretization errors which turn out 
to be large and decisive in the cases considered below. 



B. Finite size effects and anisotropics 

The pitch length of the helices is determined by the 
ratio K/J. We choose K/J = tan(27r/I0) « 0.727 
to obtain a pitch length of 10 lattice sites for a helix 
propagating in [100] direction at zero fielcP" 3 -'. We have 
found that the maximal lattice size tractable in reason- 
able CPU time is given by TV = 30 3 spins, which already 
hosts up to nine skyrmion tubes in total (However, for 
isolated cases we have checked our results against simu- 
lations on lattices of size N = 40 3 with agreeing results) . 
We use periodic boundary conditions since open bound- 
ary conditions lead to polarized spins on the boundaries 
due to missing next neighbor FM and DM interaction 
which makes them profit maximally from the Zeeman 
energy. Since it is impossible to choose parameters such 
that helices e.g. in [111] and [100] direction fit perfectly 
on the lattice at the same time one would expect strong 
finite size effects. However, we found that this turned 
out not to be a major complication in our simulations. 

The discretization of the continuum model, on the 
other hand, creates anisotropies which have to be taken 
seriously This can be seen as follows: On the lattice, the 
FM Hcisenberg term in Eq. (pi) after Fourier transform 
reads 



H FM = J^> k S(k)-S(-k), 



(3) 



where 
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which implies that all kinds of high orders in momentum 
terms are generated (the constant term only shifts the 
energies) . If we contrast this from the Fourier transform 
of the first term in the comtinuum model, Eq. (fTl), we 
see that there only the quadratic terms are present. One 
would not be worried about the higher order terms in the 
series if the ordered state was described by a uniform spin 
texture: For instance, the critical properties of the purely 
ferromagnetic Heisenberg lattice model and the Landau- 
Ginzburg functionals are in perfect agreement with each 
other and the anisotropics do not play a role at all. This 
in general is not true, especially if the critical modes do 
not become soft at zero momentum, as is the case for 
the simple Heisenberg ferromagnet, but at a finite or- 
dering wave vector, henceforth called Q. Since we use 
relatively small lattice sizes, |Q|a in general is on the or- 
der < 1 which is not a small number. Consequently, the 
contribution of the higher order terms is not negligible 
and spoils our analysis. In order to compensate for these 
induced anisotropies we add next-nearest neighbor inter- 
actions H' to our Hamiltonian. These terms are chosen 
such that they do not break symmetries of the underly- 
ing system and give a better approximation to the con- 
tinuum field theory in the sense of rendering corrections 
from higher orders of the expansion small. They assume 
the form 

H = J 2_^ S r ■ (S r+2 x + S r+2 y + S r+ 2z) + 
r 

K' J2 ( S r X S r+2x ■ X + S r X S r+2y ■ y+ 

r 
S r X S r+2 z ' 2) . (5) 

The full ak of the Heisenberg term, see Eq. pi, is now 
given by 
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This immediately shows that we can compensate the 
anisotropies to leading order by choosing J' = J/16. Re- 
peating the same procedure for the DM term leads to 
K' = K/S. 

Another way to think about this compensation is that 
the approximation of the gradient terms in the continuum 
model, Eq. (fTl) , solely by next neighbor interactions as in 
Eq. |2]) is not accurate if the spin configuration varies 
significantly from site to site. If we could simulate larger 
lattices we could use a smaller value of K which in turn 
increases the modulation period of the helices. The spin 



configuration would then vary more smoothly and con- 
sequently the approximation in Eq. (|2| becomes better. 
To summarize, the purpose of the next-nearest neighbor 
interaction terms J' and K' is to improve the approxi- 
mation of the gradient terms in Eq. (nl) by compensating 
the relatively short pitches in our simulation. 



C. Determination of thermodynamic phases and 
MC algorithm 

The different phases in our problem can be distin- 
guished either from the real space spin textures or, more 
easily, from the spin structure factor in reciprocal space. 
We calculate the average spin configuration ( S r ) from 
usually 2000 spin configurations separated by 30 lattice 
sweeps and then Fourier transform the average configu- 
ration into momentum space, 



N ^ 



S r ) exp(— ik ■ r) 



(7) 



Afterwards, we analyze the Bragg intensity profile 
/(k) oc ||( Sk )|| 2 , which corresponds to what is mea- 
sured in neutron scattering experiments. A single helix 
with wave vector Q is characterized by two Bragg spots 
sitting at Q and — Q (as required by the real order pa- 
rameter). The helical and conical phase can thus easily 
be distinguished by the direction of Q (while Q is parallel 
to the magnetic field in the conical phase it is along [111] 
in the helical). The skyrmion lattice phase has a richer 
structure and is easily identified by its six Bragg spots 
which are arranged on a regular hexagon in the plane 
perpendicular to B. The real space spin structures for 
the skyrmion lattice phase as well as the Bragg intensity 
patterns for all three phases are shown in Fig. [I] 

We determine the spin configurations using a single- 
site Metropolis algorithm. The model at hand turns out 
to be very hysteretic, in fact we checked that even par- 
allel tempering MC (PTMC) is not able to describe the 
phase transition from the skyrmion to the conical phase. 
In order to overcome this and ensure consistency we use 
three different schemes in parallel: (i) Simulated anneal- 
ing, meaning cooling at constant field, (ii) Simulated an- 
nealing to the target temperature at zero field followed by 
slowly increasing the field, (iii) Simulated annealing to a 
target temperature at high field (such that we always re- 
main in the spin polarized phase) followed by decreasing 
the field. 

If no unique state is reached this means the single-site 
Metropolis algorithm is trapped in a metastable state. 
We then use a global update which allows the system to 
fluctuate between these states and allows to determine 
the true thermodynamic state: In practice, we simulate 
two systems in the respective different states in parallel 
and offer each system 2 • 10 5 times to take the state of the 
other one after performing 20 lattice sweeps using single- 
site updates in its own state. The transition probabilities 



p(i -► j) = min (1, exp((Ej " E t )/T)) 



(8) 



a) 

easy axis 



b) 



B 






{■ 



d) helical 




.mi 



# 



conical 




FIG. 1: a) Schematic sketch of the magnetization in the helical phase, b) Schematic sketch of the magnetization in the 
conical phase, c) Averaged magnetization in the skyrmion phase in two different crystal planes for 

(J, K, B, T) — (1, tan(27r/10), 0.18, 0.82). The magnetization in direction of the external field vanishes along the yellow tubes. 
d) Bragg intensity patterns projected into the (001) plane (which is _L B) (left) and (010) plane (right). Parameters are 
(J,K,B,T) = (l,tan(27r/10),0,0) (helical phase), (l,tan(27r/10), 0.18, 0.82) (skyrmion phase) and (l,tan(27r/10), 0.18,0) 
(conical phase). 



for these hypothetical steps are recorded and averaged. 
The detailed balance condition reads 



P(i) 

P{j) 



p(j -> i) 



(9) 



p(i -> j) 

where p(i) and p(j) are the probabilities to find the sys- 
tem in the respective states. According to RefP^, this 
can be used to calculate the free energy difference be- 
tween the two states according to 



A Fij = -Tln^- 
v{j) 



(10) 



which allows to determine the thermodynamic state by 
selecting the one with the lower free energy. 



III. GLOBAL PHASE DIAGRAM 

Our main result is the non-perturbative determina- 
tion of the phase diagram associated with the free energy 
functional introduced in Eq. |I]) , see Fig. [5] 

As mentioned in Sec. |IIB[ we have to account for the 
presence of discretization errors by subtracting the quar- 
tic terms of the nearest neighbor interaction. To show 
that our simulations are severely hampered by these ef- 
fects, we have determined the [B, T)-phase diagram with 
and without anisotropy compensation. 

In both cases our simulation finds all three ordered 
phases found in the experiment on MnSi (Fig. n] for 
their real and reciprocal space signatures). In particu- 
lar we find a helical phase for small magnetic fields even 
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FIG. 2: Phase diagram with anisotropy compensation for 
(J,K) = (l,tan(2-7r/10)). Ne xt-ne arest neighbor interactions 

as J' — J/16 and 



II B 



are chosen according to Sec. 

A" — K/8. The inset shows the phase diagram without 

anisotropy compensation. 



though our Hamiltonian does not provide any explicit 
anisotropics that favor a certain crystal direction. These 
anisotropies are, however, automatically generated in the 
lattice model due to discretization errors and seem to fa- 
vor propagation in [111] direction in our case, as explicit 
in Fig. |T] 

Most importantly, our simulation shows a stable 
skyrmion phase at intermediate fields. Fig. [2] shows the 



phase diagram obtained with anisotropy compensation. 
The skyrmion phase is stable only in a small pocket 
close to T c (c. f. Fig. [2]) as it is also observed experi- 
mentally. Without the anisotropy compensation on the 
other hand, the skyrmion phase remains stable even for 
T —¥ (c. f. inset in Fig. [2]), which is in disagreement 
with experiment. We conclude that the discretization 
anisotropics indeed spoil our analysis and the true phase 
diagram is only obtained after compensation of these ef- 
fects to leading order. The real-space spin configuration 
of the skyrmion phase obtained from our MC simulations 
is shown in Fig. [T] 

This behavior has to be contrasted from previous anal- 
ysis in the case of two dimensional thin film systems. 
In numerical simulations for two dimensional systems or 
thin films with the field perpendicular to the plane one 
did not encounter the need for anisotropy compensation. 
This is related to the fact that there the conical phase 
ceases to exist (since in the conical phase the spin texture 
likes to propagate along the magnetic field) and no com- 
petition between the conical and skyrmion phase takes 
place. Consequently, the skyrmion phase remains stable 
for T -> CP. 

To summarize, we reproduce the full phase diagram of 
three dimensional helical magnets in a non-perturbative 
manner which holds beyond mean-field plus low order 
fluctuation analysis. Our analysis conclusively shows 
that the original claim that thermal fluctuations lower 
the free energy of the skyrmion phase as compared to 
the conical phase within a finite pocket which was based 
on the lowest non-trivial order in an expansion around 
mean-fieldi'is correct and higher order corrections do not 
spoil the picture. 



IV. THERMODYNAMICS ACROSS THE 
TEMPERATURE DRIVEN PHASE TRANSITION 




FIG. 3: Specific heat and longitudinal susceptibility from 



Eq. (Ill for different magnetic fields and with anisotropy 
compensation according to Sec |IIB| Errorbars have been 
calculated using the jackknife method. An offset of 0.5 has 
been introduced to seperate different curves for the specific 
heat and To was used as a fit parameter. The experimental 
data was taken from 3 



Besides the phase diagram we have studied the temper- 
ature driven phase transition into either the helical, con- 
ical, or skyrmion lattice phase (depending on the mag- 
netic field) using PTMC. 

The thermal transition at higher fields looks like a stan- 
dard second order phase transition. At lower fields upon 
decreasing temperature there is an incipient behavior 
reminiscent of a second order phase transition. This be- 
havior is controlled by the standard Wilson-Fisher fixed 
point. Upon decreasing temperature further below the 
Wilson-Fisher scale there is a second regime in which the 
system realizes that the critical modes do not go soft at a 
point in momentum space but instead on a whole sphere. 
Consequently, there is an abundance of soft modes which 
eventually drives the transition first order. The scenario 
outlined has been put forth by BrazovskiP and was re- 
cently studied in great detail in the context of MnSi 3 . 
The thermodynamic properties of the transition obtained 
from MC show striking similarity with the experimental 
findingalSHSSl anc j the Brazovskii-scenario. We have stud- 



ied two thermodynamic quantities, the specific heat cy 
and longitudinal susceptibility \zz (B = Bz). In MC 
both these quantities can be calculated from simple av- 
erages according to 



c v (T) = 



(E 2 ) - {Ef 
NT 2 



Xzz{T) = 



(M 2 ) - {M z )< 

NT 



(11) 

Fig. [3] shows the specific heat calculated from our sim- 
ulations as well as a comparison of Xzz to experimental 
data at zero field taken from [3]. The specific heat for 
low fields clearly shows a first order peak on the low tem- 
perature side of the seeming second order transition. For 
higher fields there is a tendency of a vanishing first or- 
der peak which is an indication that the transition turns 
second order. As mentioned before, PTMC is not able to 
resolve the transition from the skyrmion to the conical 
phase and thus we do not observe any sign of this phase 
transition in our simulations. The longitudinal suscep- 



tibility compares well to the experimental data, in fact 
we find the characteristic drop of the susceptibility at Tq 
as well as the characteristic inflection point at slightly 
higher temperatures^. In all cases the MC data com- 
pares favorably to experiments. 



V. CONCLUSION 

In this paper we have analyzed a full 3 dimensional 
classical spin model that describes chiral magnets like 
e. g. MnSi on a microscopic level. This effective model 
neglects the fact that these systems in general are no 
insulators but metallic in character. We used classical 
MC to determine the phase diagram and studied thermo- 
dynamic quantities across the thermal transition, such 
as the specific heat and the longitudinal susceptibility. 



From a simulation point of view, we identified the cru- 
cial role of lattice discretization anisotropies. After ap- 
propriate compensation we were able to reproduce the 
experimental phase diagram of MnSi qualitatively and 
have conclusively demonstrated that thermal fluctuations 
alone are sufficient to stabilize the skyrmion and that 
this assertion holds beyond Gaussian level. The calcu- 
lated specific heat and longitudinal susceptibility shows 
remarkable agreement with the experimental data as well 
as recent analytical approaches. 
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